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ABSTRACT 



Numerical solutions are obtained for two-dimensional, unsteady 
potential flow about a flat plate. Two types of flow are investigated: 
small amplitude oscillation about a spanwise axis and impulsive accelera- 
tion of a nonoscillating plate. Pressure, forces and moments, and kine- 
matics of the vortex wake are calculated. 

In this work vortex sheet shedding from the trailing edge is 
included. In future efforts, vortex sheet shedding from trailing and 
leading edges are contemplated. 
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NOMENCLATURE 
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- chord length 

- drag per unit span 

- unit vectors in (x- , y- , z-) directions respectively 

- lift per unit span 

- number of time increments taken 

- number of evenly distributed bound vortices 

- position vector 

- arbitrary running variable 

- time 

- trailing vortex sheet 

- free stream velocity 

- final free stream velocity 

- field velocity at r Q 

- vortex sheet 

- absolute reference system 
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vortex strength distribution 
free vortex strength 
amplitude of oscillation 
body fixed reference system 

unit vectors in (5-, T]- , C-) directions respectively 
angular frequency 
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Nondimensional quantities 

A 
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A £ 
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Subscripts 

(x, y) - component in (x- , y-) direction 
(§ , T]) - component in (§-, T]-) direction 



/ 



2 



I . INTRODUCTION 

Present generation computers make numerical approach to unsteady, 
two-dimensional subsonic flow attractive. In this study, two types of 
flow about a flat plate in nonlinear motion are investigated. Although 
having different objectives, related studies have been carried out by 

J. P. Giesing [l], N. D. Ham [2], and F. H. Abernathy and R. E. Kronauer [3]. 

Observations show that unsteady flow about a lifting body is 
accompanied by a vortex sheet being shed from the trailing edge. In the 
present study of unsteady flow about a flat plate, the plate is also re- 
placed by a bound vortex sheet. Both the bound and the free vortex sheet 
are then replaced by a finite number of discrete point vortices; thus, 
the integral equation of potential theory representing the field velocity 
is reduced to a finite difference equation. With appropriate application, 
this finite difference equation is used to determine the positions of the 
shed free vortices as a function of time and to express the boundary 
condition at any instant of time. In regard to the boundary condition, 
it is stipulated that the strengths of the vortices, both bound and shed, 
are such that the trailing streamline is tangent to the trailing edge. 

This satisfies the Kutta condition. Together with Kelvin's theorem, 
it is sufficient to determine uniquely the strengths of all vortices at 
any given instant of time. Steady flow is treated as a degenerate case 



of the above. 
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The kinematics of the shed vortex sheet and strength distri- 
bution of the airfoil and free vortices are determined. Calculation of 
such important physical parameters as stability derivatives follows 
immediately. 

II. THE BASIC EQUATIONS 

A potential flow model is used for unsteady flow about an 
infinite flat plate. Laplace's equation governs the velocity field 
and pressure can be calculated from Bernoulli's equation. From Kelvin's 
theorem the change of circulation about the flat plate is balanced by 
the shed vortex sheet(s)q. Finally, it is assumed that the Kutta condi- 
tion is satisfied; this condition is sufficient to determine the rate 
of vortex sheet shedding. 

Using thin airfoil theory, the flat plate is replaced by a 
vortex sheet where vortex strength distribution is adjusted to make the 
flat plate a streamline of the flow. 

Consider an arbitrary flow field containing a vortex sheet 

A A 

with vortex strength distribution Y(s). Vorticity is considered positive 

A 

when inducing counterclockwise circulation, r is the position vector 



1. Bluff body flow is stipulated to be accompanied by vortex sheets 
shed from the leading and trailing edges. 
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on the sheet while s is the running variable along the sheet. 

A 

velocity at r is then given by 



The field 
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The same equation written nondimens ionally is 
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Hereafter we will omit specific reference to the temporal argument. 

A description of the geometry of the field is appropriate. 

The absolute reference system is assigned coordinates (x,y,z) with their 
origin at the geometrical center of the flat plate. In the present study 
the center of the flat plate is therefore fixed but the plate is free 
to rotate in an arbitrary manner. i, i, and k are unit vectors in the 
x, y, and z directions respectively. A body fixed coordinate system 
(5 , T|, 0 with unit vectors e^ 9 e^, and also has its origin at the 
geometrical center of the flat plate. 
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Obviously, any point in the two-dimensional space considered 
may be described by either coordinate system. Suppose £ is the vector 
from the origin to the point of interest. Then 



r «■ i X + i y 

ry* <v J * 

rst 

or 



r = e £ + e n n 

Furthermore, it is easily seen 

X = COS 0 £ - sin 0 n 
y = sin © £ +• cos 0 



Conversely, 

£ = cos © X +• sin © y 
^ - sine x + cos© y 
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If we decompose the integral into the plate portion of bound vorticity 
and the free vortex sheet portion, we can write instead of (2) 
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~ L ( y c - sin© g ) + j ( x o - Cose g) 
( X a - cosGg ) 2 + ( y o - Sin© g) 2 
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A useful result obtained from (7) is the velocity normal to 
the plate at some point on the plate (T) o = 0— ). Thus 



- e i • V(5.,o i) - 

- 1/2 

T(s) J cos© ( x o - x ) + sin© ( y 0 - y) 
2flr 1 (x o - X) 2 + (y o - y) 2 

whereas the tangential velocity is 
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ds 
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Now, for a flat plate whose angular position is given as a function of 
time, 0 (t ) , the nondimensional boundary condition to be satisfied is 




i 1 

--<£ <-+€ 
o O 0 



( 10 ) 



9 (T ) is an arbitrary input to the calculation. By extending the boundary 
condition an indefinitely small distance downstream of the trailing edge, 
the Kutta condition is satisfied; i.e„, the trailing streamline is tangent 
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to the trailing edge. To determine the strength of shed vorticity, 



one makes use of Kelvin 1 s theorem where 
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The limit of integration can include any portion of the free vortex 
sheet in a Lagrangian sense or it can include the airfoil and any 
contiguous portion of the vortex sheet. 

III. THE FINITE DIFFERENCE SYSTEM OF EQUATIONS 

The numerical scheme to arrive at solutions for the airfoil 
and free vortex strength distribution and the kinematics of the free 
vortex sheet will now be discussed. As a reference, the flow chart (see Ap- 
pendix) will be instructive. The scheme is basically two-part. Knowing 
the strengths and positions of all vortices at t = t^, the velocity 
field is calculated and the new location of the free vortices are calcu- 
lated at t = t^ + At. The circulation strength of each free vortex is 

conserved. Then the unknown strengths of the bound vortices and most 
recently shed free vortex are calculated by satisfying boundary conditions 
as described by the finite difference counterpart of (8) and (10). 



by a finite number, N, of evenly distributed bound vortices. Bound 
vortex positions are indicated by circles in Fig. 1 (where, as an example, 
we choose N = 5). The distance between the vortices is 



The vortex sheet representing the airfoil will be approximated 
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The matching points, at which the boundary conditions, equation (10), 

are satisfied will be halfway between consecutive vortices, so that as 
N -> 00 , the normal velocity will approach the principal value of the integral 
given by (8). Thus , the finite number of evenly distributed bound vortices 
will be placed at 



j - 1 

= - 05 + i (j- 1,2 N) 

‘ N - 1 J 

while the matching condition points will be located at 



(13) 



5 • - - 0 5 + l -^’ 5 U - 1.2 , N - 1 ) (14) 

ol N - 1 

The shed vortex sheet will be approximated as discrete doubly- inf inite 
free point vortices shed from the trailing edge. A new vortex is shed 
at each increment in time and its strength is uniquely determined in 
satisfying the boundary conditions (equation 10). Once the strength of 
each vortex is determined, it will remain constant for all time. 

To start the calculation, let us suppose that the vortex 

th 

sheet distribution is known. The rate of change of position of the i— 
free vortex during a time increment At is given by the finite difference 
counterpart of (7). Replacing the vortex sheet strength yAs by at 

x,y = Xj , yj 
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Then the new position of each free vortex may be determined 
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tfe have found it convenient to index by unity each free 
vortex. Thus the newly shed free vortex will always be identified by 
i = 1. Having established the positions of the free vortices according 
to (15) for t = t^. + At, combining (8) and (10) the boundary condition 
to be satisfied at that instant of time may be written in a finite 
difference sense as 
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with the right-hand side completely specified. (V^ -*- s the component 

til 

of velocity normal to the surface at induced by the j— free vortex. 

It is given by 
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It may be shown that Kelvin's theorem (11) may be written as 
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Equations (17), (18), and (19) may be combined and written as 

N*1 



Z a - t>> 

,‘J J 

J Wi 



C . (i-1, 2 ,N + l) 



l 



where 



( 20 ) 



M 



C. 

L 



l <V,I. 



- U sin© - S- - 

00 04 dt 



( 21 ) 



N*i 



and , 



N i 

I r a £ } 

J -1 t-t 



( 22 ) 



T. 

J 



(j - 1 2, N) 



(23) 



11 - 



t> 

N ♦ 



1 



r 

i 



(24) 



The coefficient matrix follows readily 
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Finally from (18) 
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Numerical representation of the boundary condition and Kelvin's theorem 
is now complete. Solution is obtained by inversion of the coefficient 
matrix. The equations were programmed in Fortran TV for the IBM 360-67 
(see Appendix). 
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IV . CALCULATED RESULTS 

The solution for the field about a flat plate in steady uniform 
flow is first considered. A comparison is made in Fig. (2) between the 
vortex strength distribution determined by the finite difference technique 
and that predicted by the exact integral (see, for example, reference [4]). 
Greater agreement is achieved by increasing N, the number of bound vortices 
defining the flat plate. This is, as noted previously, due to the fact 
that as N -> 00 , the solution from the finite difference technique approaches 
the principal value of (8). 

The case of impulsive acceleration from rest was considered 
with 0 =0.1 rad (5.73°). The resultant shed vortex sheet is depicted 
in Fig. (3) for T = 2.0 and T = 3.0. The vortex strength distribution on 
the plate approaches that of steady state asymptotically with time while 
the strength of successive shed vortices becomes negligible. Fig. (4) 
shows ci (normalized on angle of attack) vs. T. As expected, it approaches 
the analytical value of 2 tt asymptotically. 

Harmonic oscillation in pitch was then considered. Q(t) was 
defined such that 



0 (t) = 0 Q sin wt 

The important physical parameters of 0 Q and uo were set equal to 0.1 and 
0.1, 1.0, 2.0, and 3.0 respectively. The number of bound vortices 
defining the plate, N, was set to 11 and the maximum number of time steps 
taken was 150. At was fixed so that the longitudinal spacing of shed 
vortices would be of the order of A§ . 
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Figs. (5-10) show the resultant shed vortex configurations 
for the specified frequencies at selected elapsed times. Since the wake 
is approximated by discrete free vortices, curves drawn sequentially 
through the points yield an approximation to the shed vortex sheet. 

As distance downstream increases, definition of the sheet becomes more 
difficult . 

The short period frequency of a stable aircraft at cruise 
speed corresponds to a reduced frequency of the order of 0.1. As seen 
in Fig. (5), the resultant wake has negligible vertical deflection. 

This reaffirms the assumption that the wake is linear for most practical 
purposes in unsteady aerodynamics. Comparison of Figs. (6-10) demonstrates 
the increasing nonlinearity of the wake with increasing frequency. 

An interpretation is made in Fig. (10) of the vortex sheet 
for u) = 3 .0 and T = 10.0. At a sufficient distance downstream, the wake 
rolls up periodically into concentrations of free vortices. The spacing 
ratio (that is, the ratio of the vertical displacement of the "centers of 
gravity" of cloisters of vortices of opposite sign to the horizontal 
distance between cloisters of like sign) is typically 0,07. The classical 
ratio predicted by von Karraan [5] for a stable configuration of double 
infinite rows of point vortices of alternating sign is 0.281. 

Coefficients of lift and drag were calculated by applying the 
Kutta- Joukowsky Law locally to each bound vortex defining the flat plate. 
Figs. (11, 12) present these values for ^ - 3.0 as a function of elapsed 
time T. A typical value for ci at T = 0.50 and 9 = -5.71° is .28. 
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The steady state value for this particular angle of attack is .61. 
ci and are periodic with time with negligible phase shift relative 
to the harmonic oscillation. 

V. CONCLUSION AND DISCUSSION 

The method developed herein yields good results for the unsteady, 
two-dimensional flow about a flat plate executing two separate nonlinear 
motions: impulsive acceleration and small amplitude oscillation about a 
spanwise axis. With the flow field specified, such significant quantities 
as stability derivatives are obtained readily. 

In future efforts, bluff body flow, such as large amplitude 
oscillation or tumbling, with vortex sheet shedding from the trailing and 
leading edges, is contemplated. Preliminary investigations were made. 

Once again the potential flow integral equation representing the field 
velocity (1) is replaced by a finite difference counterpart. The boundary 
condition, (8) and (10), is extended an indefinitely small distance upstream 
of the leading edge and downstream of the trailing edge. (15, 17, 18, 19, 20) 
are modified accordingly. Kelvin’s theorem determines the rate of vortex 
shedding and the two Kutta conditions establish uniqueness. Initial attempts 
to establish the flow field were unsuccessful. These attempts suggested 
that further refinement is needed in the numerical model. Consideration 
must be given vortices in close proximity (possibly in the form of the 
criterion proposed by Ham [2]). Furthermore, some modification may be 
required for the boundary condition at the extremities of the lifting surface. 
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Figure 1. Coordinate Systems 
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Steady State Numerical Solution Compared to Classical 
Flat Plate Theory. 
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Figure 3. Vortex Sheet Configuration for Impulsively Started Flat Plate at Tvo Times. 




Figure 4. Coefficient of Lift (Normalized on Angle of Attack) as a Function of Time for the 
Impulsively Started Flat Plate. 




Figure 5. Vortex Sheet Configuration for Harmonically Oscillating Flat Plate. 



06 *2 



e 



00*2 06'! on l 




06*0 00*0 060 - 00 * 1 - 
k 



B 

d 

a 



B 

w 

8 

*d 

B 

d 

8 

B 

V-* 

8 



B 

d 



8 

ui 



B 

«.» 



fix C\J 
ID 

cn 

OJ 

B 

d 

in 



8 



B 

d 



e 

•n 



II 

UJ 

CJ 



CD 

C_J 



B 

oi 



O 

o 



8 

oi 



M 



B 

8 



O 



UJ 

ZD 

C3 

UJ 

oc 



CD 

UJ 

B CJ 

d ID 

CD 
UJ 
OC 

8 



O 

B ° 
O 
II 

8 

♦ Ij, k— 

os*t- no'?- 



Figure 6. Vortex Sheet Configuration for Harmonically Oscillating Flat Plate 




Figure 7. Vortex Sheet Configuration for Harmonically Oscillati 
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Figure 10. Vortex Sheet Configuration for Harmonically Oscillating Flat Plate 
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Figure 11. Coefficient of Lift as a Function of Time for Harmonically Oscillating 
Flat Plate With U) = 3.0. 
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Figure 12. Coefficient of Drag as a Function of Time for Harmonically Oscillating Flat 
Plate With <*> = 3.0. 
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OSCILLATING FLAT PLATE 



COMPUTER PROGRAM 



ooo oooo ooooooooooo ooo 



PUNCH DcCK, OSCILLATING PLAT PLATE, 



RE AL* 4NXM , NYM 

DIMENSION E(Il,IL), F ( 1 1 I , NG(ll), NP(ll), Hill), (J ( l 1 ) 

DIMENSION XB (11), YB(II), XH200), VH200), XM ID , YM(ll) 
DIMENSION uAMMAF (200), XLM(Il) 

DIMENSION A { 12, 12), 6112), CU2), L(I2), NV(12), VNM(IC) 

DIMENSION T GAM (ZOO) , 2XI2G0), ZY(2G0), X 1 B ( 1 1 ) 

PI = 6.28316 

N IS THt NUMBER OP EVENLY DISTRIBUTED BOUND VORTICES LN 1 HE LIPP- 
ING SURFACE. AN IS THE INCREMENTAL DISIANCE IN IHE PIMIE DIFFE- 
RENCE SCHEME. 

T IS THE ELAPSED TIME, TT IS THE INITIAL TIME AND DcLT IS THE TIME 
INCREMENT . 

IHETAZ IS THE AMPLITUDE OP OSCILLATION IN RADIANS WHILE FREE IS 
THE REDUCED FREQUENCY OP OSCILLATION. 

N= 1 1 

AN = 1 . 0/ FLuAT ( N— 1 ) 

N1 = N-i 
N2 = N+l 
THETAZ = 0.1 
FREQ = 0.1 
CELT = C.l 
T = 0.0 
TT = -1.0 
TZ = T - TT 
DU 130 J = 1, N 

X 1 8 ( J ) = -0.5 + AN * PLOAT(J-l) 

X 1 M { J ) = -0.5 + AN * (0.5 + FLOAT ( J- 1 ) ) 

130 CONTINUE 

THETA IS IHE MODIFIED FREE STREAM ANGLE WHILE C1HETA IS THE CHANGE 
IN THETA WITH RESPECT TO TIME. 

ThETA = THETAZ * SIN(PREG * TZ) 

DEGREE = THETA * 57.29583 

DTHETA = IHtTAZ + FRlQ * CCS ( FREQ * TZ) 

COSTii = CO S ( I HETA ) 

SINTH = S IN ( ThETA) 

XTE = 0.5 * COST H 
YTE = 0.5 * SINTH 



XXX = 


SCRT( XTE 


£ 


XTE + YTE * YTE) 


NXM = 


-YTL/XXX 






NYM = 


XTE/XXX 






Du 14 0 


J = 1 , N 






Xo ( J ) 


= X IB ( J ) 


*• 


COS TH 


Y b ( J ) 


= X 1 6 ( J ) 


* 


SINTH 


X M ( J ) 


= X1M(J) 


* 


COSTH 


YM( J) 


= X 1M ( J ) 




SINTH 


CONTINUE 






DO 150 


I = 1 , N 






DO 15 0 


J = 1 , N 






E ( I , J ) 


= 1 . / ( P I 




* (FLOAT ( J- 1 ) - 0 



150 CONTINUE 



o o o o oooooo oooooooo ooo oooooo 



c 

c 

c 

c 

c 

Du 35 o i = 1 1 N 
i- ( i ) = NXM 
j>50 CONTINUE 
1 R U W — 1 

CALL M I N V ( E t N t D » N G » N P ) 

CALL GMPRD(E,F,HiN,N,IRGW) 

DU GO I - l » N 
3(1) = h( 1 ) 

60 CONTINUE 

Iht INITIAL BOUND VORTEX STRENGTH DISTRIBUTION IS SPcCIFICD. 
CGNTINUc 

ThE UNSTtADY FINITE DIFFERENCE SGLUTICN FOLLOWS. 

DO 1GC0 M= 1 » 150 
T = T + DELI 
TZ = T - TT 

ThE FREE VORTICES ARL ASSUMED SHED FRGM THE TRAILING cDGE. 

X F ( 1 ) = X B ( N ) 

YF( 1) = Yb(N) 

ThE VELOCITY FIELD IS DETERMINED. 

CD 800 J = 1, M 

VELOCITY INDUCED BY THE I TH FREL VORTEX GRIGINALLY SHED FRLM ThE 
TRAILING EDGE GN THE J TH FREE VCRTEX GRIGINALLY SnED FKCM THE 
TRAILING EDGE. 

VXTT = C.G 
VYTT = 0.0 
IFlM.EG.l) oO TO 501 

GAMMA F ( 1 ) IS AS YeT UNKNOWN . 

DO 50 0 I = 2 » M 

A VORTEX CANNOT ACT GN ITSELF. 

IF( I. EC. J ) GO TO 5 GO 
X = XF ( J ) - XF ( I ) 

Y = Y F ( J ) - YF( I ) 

XX = X * X 
YY = Y * Y 

DEN CM = PI * (XX + YY) 

VXTT = ( ( — Y * GAMMAF ( I )) /DENCM) + VXTT 
VYTT = (( X * GAMMAF ( I ) ) /DENOM ) + VYTT 

500 CONTINUE 
501 CGNTINUc 
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VELOCITY INDUCED BY THE ITH BOUND VORTEX UN THE JTH FREE VORTEX 
UK I G I NALL Y SHED FROM THE IRAILING EDGE. 

VXBT = C.O 
VYBT = 0.0 

THE BOUND VuRTEX AT XB(N) IS EXCLUDED TO ARTIFICIALLY SOPPRtSS 
THE SINGULARITY. 

CO 60 0 I = 1 1 N 

IF( I.EQ.N.AND.J.EQ.l) GO TO 600 
X = X F ( J ) - XB ( I ) 

Y = Y F ( J ) - Yfe ( I ) 

XX = X * X 
YY - Y * Y 

DE.MGM = PI * (XX + YY) 

VXB I = ( ( -Y * B ( I ) * AN ) /DEN GM ) + VXBT 

VYBT = (( X * B(I) * AN ) /DENCM ) + VYBT 

oOO CONTINUE 

THE POSITIONS OF THE FREE VORTICES ARE DETERMINED FCR 
T = T ( K ) + DEL T . 

ZX(J) = ( V X T i + VXDT) * DcLT + DELT + XF(J) 

ZY ( J ) = ( VYTT + VYBT) * DELT + YF ( J ) 

SOC CONTINUE 

DU Bo 0 J = I » M 
X F ( J ) = Z X ( J ) 

Y F ( J ) = ZYIJ) 

660 CONTINUE 

THETA = THETAZ * S I N ( FREQ * TZ) 

DEGREE = THETA * 57,29663 
ALPHA = - ( DcGREt ) 

DTHETA = THETAZ * FREQ * COS(FREU * TZ) 

COSTH = COS l THE TA ) 

SINTH = SIN l THETA) 

XT E = 0.5 * COSTH 
YTE = 0.5 * SINTH 
XL E = — X T E 
YLE = -YTE 

XXX = S (. R T ( X T E * XTE + YTE * YTE) 

NXH = -YTE/XXX 

NYM = XIE/XXX 

Du 230 J = 1 , N 

X 6 ( J ) = X 13 ( J ) * COSTH 

Y 6 ( J ) = XiB(J) * SINTH 

X M ( J ) = XLM(J) * COSTH 

YM(J) = XIM(J) * SINTH 

230 CONTINUE 
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CALCULATING The COMPONENTS UF THE ROW VECTOR, C ( 1 ) 



C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 

C EACH ELEMENT OF THE PC* VECTOR IS, AT A PARTICULAR MAICHlfw 

C CONDITION POINT, AN EXPRESSION OF TnE NORMAL COMPONENTS CF FREE 

C STREAM VELOCITY AND VELOCITY INDUCED BY FREE VCKT1CIES (EXCEPT 

C THOSE MOST RECENTLY SHED) PLUS I HAT CGMPCNENT CUE IC ROTATION CF 

C TFc LIFTING SURFACE. 

C 

DC 200 I = I , N 

CALCULATING THE INDUCED VcLOCITY AT THE I TH MATCHING CONDITION 
POINT DOE TO THE J TH FREE VORTEX ORIGINALLY SHED FROM THE TRAILING 
EDGE, WITH ThE EXCEPTION OF THAT ONE MOST RECENTLY SHEC. 

VX = C.O 
VY = C.O 

IF (M.Eu.l) GO TO 101 
DO 100 J = 2, M 
X = XM( I) - XF ( J ) 

Y ^ Y M ( I) - YF( J ) 

XX = X * X 
YY = Y * Y 

DENGM = PI * (XX + YY) 

VX = ( ( — Y * GAMMAF ( J ) ) /DENGM ) + VX 
V Y - ( ( X * GAMMAF ( J )) /DENQM) + VY 
100 CONTINUE 
101 CONTINUE 

VNM(I) = (VX * NXM + VY * NYM) 

C (I) = NXM + VNM(I) - (X1M(1) * DTHETA) 

200 CONTINUE 

COMPONENT DUL TO KELVIN'S THEORtM. 

SUMGAM = 0.0 
DO o2G I = 1,N 
SUMGAM = 6(1) * AN + SUMGAM 
620 CONTINUE 

C ( N 2 ) = SUMGAM 

CONSTRUCTING THE COEFFICIENT MATRIX. 

ThOSE ELEMENTS DUE TO THE VELOCITY INDUCED BY The JTh BOUND VERTEX 
AT THE JTH MATCHING CONDITION POINT. 

00 210 I = 1 , N 

CO 11C J = 1, N 

A ( I , J ) = 1 . / ( P 1 * (FLOAT(J-I) - 0.5)) 

110 CONTINUE 
2 1C CONTINUE 



r 
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IHOSE ELEMENTS DUE TU THE EEEECT OF THE FREE VERTEX HOST RECENTLY 
SHED FROM THE TRAILING EDGE CN THE I TH MATCHING CCNDITICN PCINT. 

DO 93 0 I = I » N 
X = XM{ I ) - XF ( L ) 

Y = Y M { I ) - YF ( 1 ) 

XX = X * X 
YY = Y * Y 

DENCH = PI * (XX + YY) 

VX = -Y/DENGM 
V Y = X/CENOM 

A ( I , N + 1 ) = - ( N X M * VX + NYM * VY) 

530 CONTINUE 

THOSE ELEMENTS DUE TO CONS I D ERA T I ON OF KELVIN'S ThECREM. 

CO 550 J = 1,N 
A ( N 2 » J ) = AN 
55 0 CONTINUE 

A { N 2 , N 2 ) = I. 



INVERSION OF A( I , J ) . 



CALL M l N V ( A , N 2 , D » L » N V ) 

INVERSE { A ( I,J) ) * C ( I > = 3 ( J ) 

CALL GMPRD( A,C,6,N2,N2, IMGR) 

GAMMA F ( 1 ) = C( N2 ) 

IF(M.EG.lbO) GO TO 999 
GO TO 699 
999 CONTINUE 

RK I TE ( 7 ,20 ) T, FREQ, ALPHA 

20 FORMAT (2r 3. 2, F10.5) 

R R I T t ( 7,23) XLE, YLE 

25 FORMAT ( 2F 10 . 5 ) 

W R IT E ( 7 , 2 1> ) XTE , YTE 

26 FORMAT ( 2F 10. 5) 

»i A IT E ( 7 » 2 1 ) ( XF ( J ) , YF ( J ) , J = 1,M) 

21 Format ( 2F io. 5) 

WK HE ( 7,22) 

22 FORMAT ( • 1C CO . • ) 

699 COwTINUe 

DO 700 I = I , M 
TGAM(I) = GAMMAF ( I ) 

ZX ( I ) = XF ( I ) 

ZY { I ) - Y F ( I ) 

7co Continue 

DO 710 i = l , M 
GAMMA F ( I + I ) — T GAM ( I ) 

X F ( I + 1 ) = ZX ( I ) 

Y F ( I + I ) = Z Y ( I ) 

7io Continue 
1CC0 CCNTINUc 
STOP 
END 




thesW688 



A ZZu iCal me,hod 10 determine the two- 




3 2768 001 90581 3 
DUDLEY KNOX LIBRARY 


























